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Abstract. The radiation hydrodynamic code C05BOLD has been supplemented with the time-dependent treatment of chemical 
reaction networks. Advection of particle densities due to the hydrodynamic flow held is also included. The radiative transfer 
is treated frequency-independently, i.e. grey, so far. The upgraded code has been applied to two-dimensional simulations of 
carbon monoxide (CO) in the non-magnetic solar photosphere and low chromosphere. For this purpose a reaction network has 
been constructed, taking into account the reactions which are most important for the formation and dissociation of CO under 
the physical conditions of the solar atmosphere. The network has been strongly reduced to 27 reactions, involving the chemical 
species H, Flo, C, O, CO, CF1, OH, and a representative metal. The resulting CO number density is highest in the cool regions of 
the reversed granulation pattern at mid-photospheric heights and decreases strongly above. There, the CO abundance stays close 
to a value of 8.3 on the usual logarithmic abundance scale with [H] = 12 but is reduced in hot shock waves which are a ubiquitous 
phenomenon of the model atmosphere. For comparison, the corresponding equilibrium densities have been calculated, based on 
the reaction network but also under assumption of instantaneous chemical equilibrium by applying the Rybicki & Hummer (RH) 
code b v lUitenbroekl fjoPll) . Owing to the short chemical timescales, the assumption holds for a large fraction of the atmosphere, 
in particular the photosphere. In contrast, the CO number density deviates strongly from the corresponding equilibrium value in 
the vicinity of chromospheric shock waves. Simulations with altered reaction network clearly show that the formation channel 
via hydroxide (OH) is the most important one under the conditions of the solar atmosphere. 
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1. Introduction 

Since lNoves & Hallh 1 972b ) deduced very low brightness tem¬ 
peratures in the (upper) photosphere from observations in the 
core of the carbon monoxide (CO) 3-2 R14 line, the thermal 
structure of this layer has been subject of an ongoing controver¬ 
sial debate. In many subsequent observations the low tempera- 
tures h ave not only been confirmed, e.g. bv lAvres & Teste rmanl 
dl98 lh . but also found to be a s low as 4000 K and even d own to 
only 3700 K (see Figs. 4-5 in lAvres~& Testermanil98lb . These 
small values fall below the values implied by c omm only used 
semi-empirical models like VAL-C (fVernazzaet~allll98 ll ) that 
are based on other diagnostics, e.g., in the UV range. Moreover, 
the CO observations did not show a prominent temperature in¬ 
version which is an important feature of VAL-C-like models. 
Ayres & Testerman explained this ’diagnostic dilemma’ with 
the thermal bifurcation of the outer solar layers. Following that 
idea, the hot plasma would be confined in discrete structures 
embedded in cool gas that can account for the observed carbon 
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monoxide features. The inhomogeneities would be preferably 
of small spatial scale and large thermal contrast. 

A number of investigations imply that the bulk of CO is lo¬ 
cated in the low chromosphere and below. For instance, Ayres 
& Testerman state that the CO concentration peaks at an opti¬ 
cal depth of T 500 ^ 1CT 2 for a reference wavele ngth of 500 nm. 
That is in line with the results hv LSolanki et alj ( 122.4 who ob¬ 
served CO in emission at the solar limb. There, it remains at 
a constant level up to 0"4 (~ 300 km) above the (continuum) 
limb bef ore it decreases rapidly in the layers above (see their 
Fig. 2). llJitenbroek et alj (1 994 ) carried out similar observa¬ 
tions, but with an improved technique, and found CO in emis¬ 
sion, extending 075 (« 360 km) beyond the continuum limb. 
lAvres & Rabin 0996 ) derive, based on the AV = 1 bands, an 
off-limb extension of 0'.'6 (~ 440 km) for the strongest lines, 
indicating the presence of cool gas up to 350 km above the 
classical temperature minimum. 

Given the high temperatures implied by other diagnostics, 
the above-mentioned CO observations inevitably demand for 
a spatially inhomogeneous structure of the photosphere and 
low chromosphere. In this sense. Solanki et al. (1224), based 
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on observed horizontal velocities, suggest that the bulk of CO 
might be located above granule interiors. Direct evidence for 
this assumption was provided by llJitenbroeki ( 2000 ., see also 
Uitenbroek et al., 1994), who took spectroheliograms in the 
cores of CO lines. They discovered not only a pattern con¬ 
nected to the magnetic network, but also bright rings with dark 
centres in the quiet Sun, resembling a reversed granular pattern. 
Moreover, they conclude that such spatial variations are largely 
of dynamic nature, and that dynamics thus play an important 
role in the formation of the dark CO line cores. In fact, oscilla¬ 
tory intensity changes_(see, e.g. jAvresl2004 have already been 
detected bv lNoves & Halil J 1 972bl> . 

A number of static models have been constructed in or¬ 
der to explain the observations, ranging from pure ID at¬ 
mospheres like the COmos nhere JWiedemann et al.lll994l) to 
two-c omponent models (e.g. lAvres et alJll 9861 1 Avres & Rghin 
1996) and ev en more complex spatial configurations (e.g. lAvres 
199 li . 2002 1. But although indicating the right direction, the 
obvious dynamic nature of CO features cannot be reproduced 
by means of a static approach. Spatial and temporal variations 
have to be taken into account for a detailed model. 

Instead of inverting observed intensities_to_a temperature 
stratification with one or more components, llJitenbroek 1 2000 11 
started from a_ snapshot of the 3D hydrodynamic al model by 
Stein & Nordlund: (T 989I ) and calculated CO number densities 
under the assumption of instantaneous chemical equilibrium 
(ICE). The highest concentration was found in the middle pho¬ 
tosphere (~ 100 - 300 km above optical depth unity) above 
granule interiors, while the concentrations above intergranu¬ 
lar lanes (at the same height) are 3 - 4 times smaller. The CO 
distribution is strongly connected to the cooling due to strong 
adiabatic expansion and cooling above the granule centres. The 
calculations were performed under the assumption of instanta¬ 
neous chemical equilibrium. Neither the radiative cooling ac¬ 
tion by CO itself nor advection due to the hydrodynamic flow 
were taken into account. This pioneering work showed the im¬ 
portance of solar granulation for CO line formation. 

Asensio Ramos et al. (1200 31 proceed to non-equilibrium 
CO chemistry calculatio ns based on the tim e-dependent hydro- 
dynamic simulations by ICarlsson & Stein (1997). Solving the 
chemistry equations for a chemical reaction network, includ¬ 
ing the most relevant species, resulted in CO number density 
as function of height and time. From this Asensio Ramos et al. 
conclude that the radiation in CO lines close to the solar limb 
originates from heights not greater than 700 km. Moreover, 
they proved that ICE is a valid assumption for the lower layers 
but overestimates the CO number densities above the low chro¬ 
mosphere. The calculations were one-dimensional only and 
thus did not allow for an analysis of the horizontal distribution. 

This paper is the first part of a series. Here we present 
time-dependent CO chemistry as part of multi-dimensional ra¬ 
diation hydrodynamics simulations, starting with the 2D case. 
The major advance with respect to earlier investigations is the 
time-dependent treatment of a chemical network in combina¬ 
tion with advection of particle densities with the hydrodynamic 
flow in two or three spatial dimensions. The upgraded code 
C0 5 BOLD dFrevtag et all2002h will serve as a base for further 


improvements as, e.g., the back reaction of CO as a cooling 
agent. 

After the description of the methods, the chemical input 
data, and the numerical model in Sects. |2J01 and^] respectively, 
we present the results of two-dimensional simulations in Sect.|^ 
which are discussed in Sect. [ 6 ] Finally, conclusions are drawn 
in Sect.|7] 


2. Numerical method 


2.1. Formulation of the problem 


The number density «, of a chemical species within a fixed vol¬ 
ume can change in time ( t ) due to advection as expressed by 
the continuity equation: 

■ 

-^ + V • (n,v) = 0 , (1) 

where v denotes the velocity of the hydrodynamic flow. An ad¬ 
ditional source term must be taken into account, if the number 
density also changes due to chemical reactions. Such changes 
can be realised by a large variety of reactions. In the present ap¬ 
plication we restrict ourselves to two- and three-body reactions. 
The new term can then be written as 


dtij 

dt 
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where n, are the number densities of the different species. The 
first and second right-hand terms present two-body reactions 
which yield losses (negative sign) and gains (positive sign) for 
species density n, with rates & 2 ,,y and kn.ji, respectively. Three- 
body reactions are analogously accounted for by the third and 
fourth term with rates £ 377 / and k^ji„,. The resulting ordinary 
differential equation Eq.[3is of first order and has to be imposed 
for each chemical species separately. The whole problem thus 
requires the solution of a system of differential equations. 

Chemical reactions can have rates that differ by many 
orders of magnitude. Hence, not only the number densities 
of different species but also their temporal derivatives cover 
a large range, causing the system of equations to be stiff. 
Consequently, an implicit scheme should be used for the nu¬ 
merical solution of the problem. 


2.2. Time-dependent solution 

The time-dependent treatment of a chemical reaction network 
has been added as an optional separate computational step 
to the radiation hydrodynamics code CO BOLD. The general 
properties of the code are described in Fre vtag et alJ 12002) and 
IWedemever et Hi! d2004l hereafter W04). Following the gen¬ 
eral concept of operator splitting, the chemistry calculations 
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are performed after the hydrodynamics step for each computa¬ 
tional time step. The treatment of advection and chemical reac¬ 
tions is thus separated into subsequent steps. First, the number 
densities of all chemical species are advected with the hydro- 
dynamic flow field. This is done analogously to the gas den¬ 
sity within the applied Roe solver. After the subsequent viscos¬ 
ity step the chemical reactions are handled for each grid cell 
separately, starting with the calculation of the reaction rates. 
Necessary input data are the local temperature, which is avail¬ 
able as result of the hydrodynamics step, and the number densi¬ 
ties of the involved chemical species of the previous time step. 
The basic rate is given by 

k = aT 3 Je~? /T , (3) 


treatment described in Sect. E3 but with the local gas tempera¬ 
ture kept constant and without advection, until the CO number 
density reaches an equilibrium value «co,eq- Starting from an 
initial value «co,o = «co ( f = 0), the resulting temporal evolu¬ 
tion of the CO number density follows in most cases perfectly 
the function 

nco(t) = (nc o,o - »co,eq) e~ ,/,cban + nco.o , (6) 

which is the typical solution for the differential chemistry equa¬ 
tion (Eq.[2}. The derived profile nco(t) then allows to determine 
the equilibrium value 

n co,CE = «co,eq = lim «co (0 (7) 

n t —>oo 


where = T /300 K. For catalytic reactions which involve 
a representative metal also the number density n M of the metal 
enters: 

k = n M a T m P e^ /T . (4) 


Given the rates, the system of differential equations is defined 
and then solved with an implicit BDF (backward differenti¬ 
ation formula^ method. Here, the DVODE package is used 
ferown et alil989i ). The solver uses an internal computational 
time step which is adjusted automatically. The solution pro¬ 
vides the number densities of the involved species after the 
overall computational time step prescribed by the foregoing hy¬ 
drodynamics. 

The radiative transfer is treated in a subsequent step 
frequency-independently (grey) and under strict assumption of 
local thermodynamic equilibrium (LTE). See W04 for more de¬ 
tails. 

The boundary conditions are consistent with the hydrody¬ 
namics part of the code, i.e. in the presented simulations chem¬ 
ical species are advected across the lateral periodic boundaries 
and can leave the computational domain at top or bottom. Like 
for the other boundary conditions, the number densities in the 
bottom grid layer are copied into the ghost cells below but are 
scaled to the mass density p of those cells: 


dpij/p) 

dz 


o 


ghost _ bottom ghost / bottom 

/ — P IP 


(5) 


This way the total particle numbers of the involved elements 
(here H, C, O, and a representative metal M) are almost per¬ 
fectly conserved, except for a decrease of typically 5 10 4 per 
simulation hour in relative number due to mass loss at the upper 
boundary. 


2.3. Chemical equilibrium 

Next to the full time-dependent simulations, separate calcula¬ 
tions have been performed in order to derive CO equilibrium 
densities. In this case, which is referred to as CE in the follow¬ 
ing, different snapshots of the time-dependent simulation se¬ 
quence are used as initial condition with the number densities 
derived from the gas density in the same way as done for the 
start model (see Sect.|4}. For each grid cell the pure chemistry 
calculations are performed analogously to the time-dependent 


and the corresponding chemical timescale f c h em for every grid 
cell. In principle the assumption of instantaneous chemical 
equilibrium (ICE) should provide the same values, although 
that approach utilises chemical equilibrium constants instead 
of a reaction network. 

For comparison, ICE number densities are calculated with 
the spectru m synt hesis code RH (“Rybicki & Hummer”) by 
ll literihroek (120011) for a number of snapshots from our time- 
dependent 2D model in the same way as done bv lllitenhroekl 
( 2000 1. The results of these calculations are here referred to as 
UICE. 

2.4. ICE and spectrum synthesis 

The RH code (see Sect. 12.3> is also used for calculating spec¬ 
tra. We have modified this code, which by default computes the 
instantaneous chemical equilibrium densities for a user-defined 
set of molecules and their constituent atoms, to also accept the 
CO densities computed with CO 5 BOLD as input. The distri¬ 
bution over the various energy levels of the CO molecule is 
assumed to be in LTE, which has been sh own to be a sound ap¬ 
proximation JAvres & Wiedemannll989k so that the CO line 
opacities and line source functions are in LTE. We compute in¬ 
tensities in the 2142 - 2145 cm -1 range (4.662 - 4.668 pm), 
which is observable from the ground and which contains CO 
lines with widely varying properties dGoorvitchl 19941) . Due to 
the non-negligible contribution of Thomson scattering to the 
total opacity at the wavelengths studied, we have to lambda- 
iterate the angle-averaged radiation field J v to convergence. 
For those computations we use the A4 angular quadrature of 
Icarlson J 1963k which has 3 rays per octant. For the intensities 
in the vertical direction a separate formal solution was subse¬ 
quently performed using the angle-averaged radiation field J v 
from the full solution to account for the scattering contribu¬ 
tions. 

Since scattering off free electrons is an important source of 
opacity in higher layers of the solar atmosphere, it is necessary 
to use realistic electron densities. The default LTE values, ob¬ 
tained by solving the Saha-Boltzmann equations for all species, 
are definitely unrealistic in the chromosphere since there the 
degree of ionisation is largely decoupled from the local temper¬ 
ature. Instead of the default LTE values we derive electron den¬ 
sities from the FAL-C model dFontenla et a 1.1 1991 ) by means 
of interpolation of the N e (Nu) dependence. 



















4 


S. Wedemeyer-Bohm et al.: Carbon monoxide in the solar atmosphere 



Fig. 1 . Chemical reaction network. The seven chemical species 
(not including the representative metal) are connected via reac¬ 
tions which are listed in Tabled The involved chemical species 
of the reactions are printed at the arrows. 


( 200(J _ http: //homepages. vub. ac. be/~akonnov). 

lAsensio Ramos etal . l2003l) . for example, based their chemi¬ 
cal network entirely on the rates from this database. 

The list of reactions is further supplemented with three re¬ 
actions present in Table 2 hv 'Avres & kahili ( 1996 1 that are not 
found in UMIST (#7000, #7001, #7002). Note that their data 
are explicitely prepared for a gas temperature of 5000 K be¬ 
cause of the authors’ interest in the CO reformation timescales 
in a previously molecule-free gas cooling below the H equi¬ 
librium temperature. 

Some reactions are available from different sources and 
their reaction rate coefficients differ noticeably. Some of the 
rates in Konnov’s database for instance are different from those 
in UMIST. They differ not only in the absolute value of the 
rate, but also in the type of reaction (temperature dependence, 
activation barrier). Some others agree surprisingly well. These 
ambiguities demand for a closer examination of the chemical 
input data. Here we present a examination of different available 
sources for reaction rates included in our chemical network (see 
TableThe different data are listed in Tableland discussed 
below. 


In order to be consistent with the chemo-hydrodynamic cal¬ 
culations in this paper, the same chemical abundances have 
been adopted (see Sect.^J. By default, the following molecules 
are taken into account in RH: Hi, Hi + , C 2 , Ni, Oi, CH, CO, 
CN, NH, NO, OH, and H 2 O. While this set is used for the 
spectrum synthesis, we did also a calculation with a reduced 
set comprising only the molecules that are present in the chem¬ 
ical reaction network (see Sect.0. The results are discussed in 
Sect. |6] 


3. Chemical reaction network 

In the present study we take into account eight chemical 
species: H, H 2 , C, O, CO, CH, OH, and M, where M stands 
for a representative inert catalyst (also shortly referred to as 
’metal’). For the standard model (see Sect.^} the abundance of 
the representative metal was set to the one of helium. Strictly 
speaking He is not a metal but the most abundant element that 
can be chosen as representative catalytic element M. This way 
an upper limit for the influence of M is provided. The species 
are connected via 27 reactions (see Table 0. Ion-molecule 
reactions and photoreactions are excluded because, follow¬ 
ing lAsensio Ramos et ali 120031) . their influence on the total 
CO concentration is negligible for heights of < 1000 km in 
the solar atmosphere. The chemical input data are taken from 
different sources 1 which are discussed below. 

First, the UMIST database rate file rate99 (Le Teuffet al. 
200 ()l) has been reduced to reactions involving chemical species 
stated above. 

In addition, a few termolecular 2 * reactions are 
taken from the combustion database of iKonnovI 

1 Note that the here given reference BDHL72 is apparently identical 
with BDDG72 as used by AR96. 

2 termolecular: A reaction involving the simultaneous collision of 

three particles. 


Reactions 8, 48. Since Ayres & Rabin fitted the UMIST 
JMillar et al.lll99ll former version) reaction rate explicitely 
for temperatures around 5000 K, their fit is poor for gas 
temperatures very different from that. The newer UMIST 
data iLe Teuff et al.l|2OO0|) comes from the NIST database 
(http: //www. nist. gov) and an accuracy of better than 25% 
for 300 K < T < 2500 K (8) and 297 K < T < 3532 K (48) 
is stated. At high temperatures, the real rates could be an or¬ 
der of magnitude larger than the UMIST rate. On the other 
hand, Konnov’s rate for reaction #48 is very similar to the one 
of UMIST. The assumption that both are probably based on 
the same experimental data supports our decision to use the 
UMIST reaction rate. 


Reaction 14. Ayres & Rabin derived this rate from the as¬ 
sumption of detailed balance for reactions #14 and #67, where 
they calculated the CO/OH equilibrium ratio from a Saha 
Ansatz. The high exponent (3 indicates that this rate is only 
valid in a very narrow temperature range around 5000 K. The 
UMIST rate is more than one order of magnitude smaller 
around 5000 K. We t race the UMIST rate back to Mitchelll 
dl984|) a nd fur ther to IWestlevI ill98(1 ). The original work of 
IWestlevI dl 98oh gives a ~ 5 times higher rate than Mitchell and 
UMIST. Since the reason for this discrepancy is unknown, we 
use the original rate from Westley. 


Reaction 6 7. For this reaction, the source of the UMIST 
database is IPrasad & Huntress! f 1 98ol ). The reaction rate did 
not change in UMIST between the 1991 and 1999 rate files. 
Konnov’s rate is — 5 times smaller than that of UMIST around 
5000 K and has no temperature dependence (8.30 10 11 com¬ 
pared to 4.5 10~ 10 ). Westley (1980) quotes a rate of 7.39 10 11 
at 5000 K, a factor of 6 lower than UMIST. Both, Konnov’s 
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Table 1 . Chemical reaction network: Involved are the species H, tb, C, O, CO, CH, OH, a representative metal (M), and photons 
(v). The corresponding reaction rates are parameterised by the coefficients a, /?, and y, following the UMIST standard (see Eqs.0 
0. The dimension of a is cm 3 s _1 and cm°s 1 for two- and three-body reactions, respectively, is dimensionless, whereas y is 
given in units of K. T he last column gives r eference s to the origin of t he rates: UM I ST: iLe T euff et a ll J200(l KCD: Konnov’s 
combustion database konnovl2000ll . W80: IWestlevl J 1 98(il BDHL72: IBaulch et alJ (1 19721) BDDG76: IBaulch et alJ J 19761) : The 
ID numbers indicate where the reactions have been found first, regardless of the final source of the coefficients which is given 
in the last column. An ID number less than 5000 refers to UMIST, whereas numbers 5000 and above are assigned to additional 
reactions found in Konnov’s database. Reactions listed bv lAvres & Rabin . C19961) are classified with numbers of > 7000. 


ID nr. 

reaction 



a 

fi 

7 

[K] 

ref. 


radiative association 

[cm 3 s~'] 




3681 

H + C 

—» 

CH + v 

1.00(—17) 

0.00 

0.0 

UMIST 

3683 

H + O 

—» 

OH + v 

9.90(—19) 

-0.38 

0.0 

UMIST 

3707 

C + O 

-> 

CO +v 

1,58(—17) 

0.34 

1297.4 

UMIST 


3-body association 

[cm 6 s” 1 ] 




5001 

H + H + H 2 

—> 

H 2 + h 2 

9.00(-33) 

-0.60 

0.0 

KCD 

5002 

H + H + H 

—> 

H 2 + H 

4.43(-28) 

-4.00 

0.0 

BDHL72 

7000 

O + H + H 

—> 

OH + H 

1.00C-32) 

0.00 

0.0 

BDHL72 

7001 

C + O + H 

—> 

CO + H 

2.14(—29) 

-3.08 

-2114.0 

BDDG76 


Species exchange 

[cm 3 s~'] 




1 

H + CH 

—> 

c + h 2 

2.70(—11) 

0.38 

0.0 

UMIST 

8 

H + OH 

—> 

o + h 2 

6.99(—14) 

2.80 

1950.0 

UMIST 

14 

H + CO 

—> 

OH + C 

5.75C-10) 

0.50 

77755.0 

W80 

42 

h 2 + c 

—» 

CH + H 

6.64(—10) 

0.00 

11700.0 

UMIST 

48 

h 2 + o 

—> 

OH + H 

3.14(—13) 

2.70 

3150.0 

UMIST 

66 

C + OH 

—> 

O + CH 

2.25(—11) 

0.50 

14800.0 

UMIST 

67 

C + OH 

-> 

CO + H 

1.81(—11) 

0.50 

0.0 

W80 

102 

CH + O 

—> 

OH + C 

2.52C-11) 

0.00 

2381.0 

UMIST 

104 

CH + O 

—> 

CO + H 

1.02(-10) 

0.00 

914.0 

UMIST 


collisional dissociation 

[cm 3 s~*] 




4060 

h + h 2 


H + H+H 

4.67(-07) 

-1.00 

55000.0 

UMIST 

4061 

H + CH 

—> 

C + H+H 

6.00(-09) 

0.00 

40200.0 

UMIST 

4062 

H + OH 

—> 

O + H+H 

6.00(-09) 

0.00 

50900.0 

UMIST 

4069 

H 2 + h 2 

—> 

H 2 + H+ H 

1,00(—08) 

0.00 

84100.0 

UMIST 

4070 

H 2 + CH 

—> 

C + h 2 + h 

6.00(-09) 

0.00 

40200.0 

UMIST 

4071 

H 2 + OH 

—» 

O + H 2 + H 

6.00(-09) 

0.00 

50900.0 

UMIST 

7002 

CO + H 

—> 

C + O+H 

2.79(-03) 

-3.52 

128700.0 

BDDG76 


collision induced dissociation 

[cm 3 s~'] 




4076 

CO + M 

—> 

O + C+M 

2.79(—03) 

-3.52 

128700.0 

BDDG76 


catalysed termolecular reactions 

[cm 6 s" 1 ] 




4079 

H + M + O 

—> 

OH + M 

4.33(—32) 

-1.00 

0.0 

UMIST 

5000 

H + M + H 

—> 

H 2 + M 

6.43C-33) 

-1.00 

0.0 

KCD 

4097 

C + M + O 

—> 

CO + M 

2.14(—29) 

-3.08 

-2114.0 

BDDG76 


rate and the one of Westley are valid for typical combustion 
temperatures, while the UMIST rate is only valid up to 300 K. 

There are two reasons for choosing the lWestlev ( 1198 0!) rate: 

(1) the reactions #14 and #67 form a pair of forward and back¬ 
ward reaction and hence the same data should be used for them. 

(2) This reaction rate is valid for typical combustion tempera¬ 
tures and hence more appropriate to the temperature range cov¬ 
ered in our solar atmosphere simulations. 


Reaction 3683. Ayres & Rabin state UMIST as source for this 
rate coefficient although apparently they slightly altered it. For 
the present chemical network the original UMIST rate is used. 


Reaction 3707. After detailed examination of the original 
data from fDaluarno et al . ill 99(1 and the fits by Ayres & 
Rabin and UMIST, we adopted the rate coefficients provided 
by UMIST dLe Teuff et alJUOOol) . The latter give the best fit in 
the temperature range between 2000 and 8000 K. 

Reaction 4060. Again, Ayres & Rabin’s fit of the rate by 
IBaulch et al J ( 1972 ) holds around T ~ 5000 K. Outside that 
temperature range it systematically overestimates the original 
data by orders of magnitude. In contrast, the UMIST rate de¬ 
viates by not more than a factor of 3 from the original data 
in the range 3400 K to 5000 K. Even though the source of the 
UMIST fit is unknown to us, we favor it because the rate covers 
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Table 2. Comparison of chemical data from different sources. The reaction rates are parameterised by the coefficients a, /?, and 
y, following the UMIST standard. The dimension of a is cm ’s 1 and cm 6 s _1 for two- and three-body reactions, respectively. /3 is 
dimensionless, whereas y is given in units of K. The sixth column gives the temperature range in K over which the reaction rate 
is valid and the seventh column its accuracy (following the UMIST error nomenclature. A: < 25%, B: < 50%, C: within a factor 
2, D: within an orde r of magnitude, E: high ly unc ertain). The last column gives referenc es to the origin of th e rates: U MIST: 
iLe Teuff etaL dZOOOl), A R96: Avres & R abinldl996l). K CD: Ko nnov’s combustion database dKonno vl2000l) , W80: IWestlevldl 98(jl) 


BDHL72~ lBaulch et alJdl972) BDDG76: lBaulch et alJdl976t) 


ID nr. 

reaction 



a 

(3 

7 

T range 

acc. 

ref. 

1 

CH + H 

-» 

c + h 2 

2.70(—11) 

0.38 

0.0 

300 - 2000 

B 

UMIST (1) 





1.31(—10) 

0.0 

80.5 



KCD 

8 

H + OH 

-> 

o + h 2 

6.99(—14) 

2.80 

1950.0 

300 - 2500 

A 

UMIST (8) 





3.35(-13) 

1.70 

0.0 

-5000 


AR96 

14 

H + CO 

-> 

OH + C 

l.lO(-lO) 

0.50 

77700.0 

2590-41000 

C 

UMIST (14) 





3.21(—35) 

16.20 

0.0 

-5000 


AR96 





5.75(-10) 

0.50 

77755.0 

typical T comb 


W80 

48 

o + h 2 

—> 

OH + H 

3.14(—13) 

2.7 

3150.0 

297 - 3532 

A 

UMIST (48) 





2.86(—13) 

1.90 

0.0 

-5000 


AR96 





4.10(—13) 

2.7 

3165.1 



KCD 


67 

C + OH 

-» CO + H 

l.lO(-lO) 

0.5 

0.0 

10-300 C 

UMIST (67) 




1,22(—10) 

0.5 

0.0 

-5000 

AR96 




8.30(—11) 

0.0 

0.0 


KCD 




1.8K-11) 

0.5 

0.0 

typical T comb 

W80 


102 

CH + O 


C + OH 

2.52(—11) 
2.52(—11) 

0.0 

0.0 

2381.0 

2380.1 

10 - 6000 

B 

UMIST (102) 
KCD 

103/4 

CH + O 

-> 

CO + H 

6.60(—11) 

0.0 

0.0 

10 - 2000 

A 

UMIST (103) 





1,02(—10) 

0.0 

914.0 

2000 - 6000 

A 

UMIST (104) 





6.64(-ll) 

0.0 

0.0 



KCD 

3683 

H + O 

-> 

OH + hv 

9.90(—19) 

-0.38 

0.0 

10-300 

C 

UMIST (3683) 





9.24(—19) 

-0.40 

0.0 

-5000 


AR96 

3707 

C + O 

—> 

CO + hv 

1,58(—17) 

0.34 

1297.4 

300- 13900 

B 

UMIST (3707) 





5.551-18) 

0.60 

0.0 

-5000 


AR96 

4060 

H + H 2 

—» 

H + H + H 

4.67(-07) 

-1.00 

55000.0 

1833-41000 

C 

UMIST (4060) 





3.45 (-21) 

6.60 

0.0 

-5000 


AR96 





8.86(-04) 

-4.00 

51900.0 

3400 - 5000 

D 

BDHL72 

4062 

H + OH 

-> 

O + H + H 

6.001-09) 

0.00 

50900.0 

1696-41000 

C 

UMIST (4062) 





3.921-26) 

10.40 

0.0 

-5000 


AR96 

4076 

CO + M 


O + C + M 

2.861-03) 

-3.52 

112700.0 

2000 - 10000 

B 

UMIST (4076) 





2.791-03) 

-3.52 

128700.0 

7000- 15000 

D 

BDDG76 

5000 

H + H + M 

-> 

H 2 + M 

6.431-33) 

-1.0 

0.0 



KCD 





5.881-33) 

-1.0 

0.0 

1700 - 5000 

C 

BDHL72 

5001 

H + H + H 2 

— » 

H 2 + h 2 

9.001-33) 

-0.6 

0.0 



KCD 





2.391-32) 

-1.0 

0.0 

2500 - 5000 

D 

BDHL72 


5002 H + H + H 

-> H 2 + H 

8.82(-33) 

0.0 

0.0 

KCD 



4.631-28) 

-4.0 

0.0 -5000 

AR96 



1.831-31) 

-1.0 

0.0 

P83 



4.43(-28) 

-4.0 

0.0 3400 - 5000 D 

BDHL72 


7000 

O + H + H 

-> OH + H 

1.001-32) 

2.761- 33) 

2.761- 32) 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

-5000 

low.limit 
upp. limit 

E 

E 

AR96 

BDHL72 

BDHL72 

7001 

C + O + H 

-» CO + H 

3.781-29) 

2.141-29) 

-3.50 

-3.08 

0.0 

-2114.0 

-5000 

7000 - 14000 

D 

AR96 

BDDG76 

7002 

CO + H 

-> C + O + H 

1.391-46) 

2.791-03) 

22.80 

-3.52 

0.0 

128700.0 

-5000 

7000 - 15000 

D 

AR96 

BDDG76 


a larger temperature range and nicely follows the trend of the 
experimental d ata, which display a steeper gradient — steeper 
than the lBaulch et al. ( 1972 1 formula— at temperatures below 
3400 K. Baulch et al. state that their rate is only a tentative 
suggestion in a narrow temperature range. 


Reaction 4062. The Ayres & Rabin fit is again only valid for 
temperatures around 5000 K. We use the UMIST data instead, 
because it covers the whole temperature range of interest. 
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Rea ction 4076 . The UMIST database refers to 
IPetuchowski et all dl989l) . and from there to iBaulch et alJ 
Tl97~4 (see rate #7002). Apparently the exponent ft changed in 
the referencing process, the factor a suffered a rounding error 
and the error margin has been quo ted wrongly. Hence, we stick 
to the original rate of lBaulchet ai1 dl97d) . 

Reaction 5002. We cannot judge the accuracy of Konnov’s 
rate over the temper ature interval 2000 K < T < 8000 K 
but he refers back to ICohen & Westheri] dl 983l) . IPalla et aTI 
(1983) refer back to a publicatio n, where Cohen is a co-auth or 
d.Tacobs et alJll967h . In this case. ICohen & Westberd d 19831) is 
newer and presents a literature review for rate constants. The 
difference in the rates by Konnov and iBaulch et al.1 ( 1972 ) is 
the strong temperature dependence of the rate in the latter ref¬ 
erence. Since Baulch et al. (1972) is a well established and doc¬ 
umented database - often used by UMIST as well -, we use 
their rate constants. 

Reaction 7000. Owing to the lack of alternatives we adopt the 
rate hv IBaulch et al.1 ( 1972 ) altough it should be regarded as a 
mere order of magnitude estimate. 

Reactions 7001, 7002. Ayres & Rabin’s rate i s ag ain re¬ 
stricted to the 5000 K temperature range. The IBaulch et al.1 
{mi At assumes that H is very similar to Ar and it covers 
a larger temperature range, 7000 K < T < 15 000 K. Hence, 
we use the latter rate. 


4. Numerical model 

The two-dimensional numerical model consists of 120 by 
140 grid cells each with a horizontal width of 40 km. The 
height of the cells is 50 km at the bottom (at z = -1484 km) 
and smoothly decreases to 14 km for all heights above z = 
-531 km. The upper boundary is located at 1016 km. Note that 
the origin of the height axis is defined by the horizontal and 
temporal average of Rosseland optical depth unity (r = 1). The 
total extent of the computational domain is thus 4800 km x 
2500 km. Most numerical parameters are adopted from the re¬ 
cent 3D simulations by Wed emever et alJ ( |2QQ4|), including the 
usage of grey OPAL-P HOENIX opacity Jlglesias et alJll992t 
lHauschildt et al.lll997l) . The preliminary start model was ex¬ 
tracted from the 3D model, too. In order to ensure relaxation 
the model was advanced several simulation hours before it was 
finally supplemented with arrays for number densities of the 
involved chemical species. For each grid cell the same constant 
chemical composition is assumed. The required abundances 
of the a tomic species are set to fCl — 8.39 iAsnlund et all 
lin nressh . and [O] = 8.66 iAsnlund et alJl2004l) . where [H] = 
12. The initial abundances of the molecular species are set to 
1 0 211 times the hydrogen density for CO, CH, and CO, and to 
1 0 4 for Hi. The abundance of the representative metal is set 
equal to the helium abundance ([He] = 11.00) and thus pro¬ 
vides an upper limit for the influence of the metal. The initial 
number densities are then directly calculated from the gas den¬ 


sity of the final start model for the constant chemical composi¬ 
tion described above. 

For the standard model in total 86000 s have been cal¬ 
culated. The first 36000 s are reserved to ensure a sufficient 
chemical relaxation of the model and are therefore excluded 
from analysis. The global time step, which is relevant for hy¬ 
drodynamics and radiative transfer, was typically 0.2 s - 0.4 s, 
whereas the time step for the chemistry calculations has been 
adjusted by the solver itself for each individual cell and global 
time step. 

In addition to the standard model more simulations with al¬ 
tered chemical network have been produced in order to investi¬ 
gate the influence of individual reactions or reaction groups on 
the formation of CO. All simulations covered a time span of at 
least 7000 s. See Sect. I5.5l for more details. 


5. Results 


5.1. Two-dimensional distribution 

In Fig. Qa snapshot of the 2D model after a simulation time of 
70390 s is presented. Each panel displays a different quantity as 
function of horizontal ( x ) and vertical (z) position for the time- 
dependent simulation (TD) but also for the equilibrium calcu¬ 
lations CE and UICE (see Sect. l2.3> . all based on the same time 
step . The first panel displays the gas temperature, exhibiting a 
few granules and intergranular downflows. In contrast to posi¬ 
tions above granule interiors the temperature is increased above 
the downflows. These high-temperature regions in the middle 
photosphere, which form the reversed granulation patt ern when 
seen from above (see iLeenaarts & Wedemever-Bohrd 2005, 
W04), are due to compression heating as result of convective 
overshoot. There are also high-temperature regions in the up¬ 
per layers which are produced by propagating shock waves 
(see W04). An example for the latter can be seen in the up¬ 
per part of panels a and b, outlined by the contour line for 
T = 5000 K. Both phenomena induce spatial and temporal 
inhomogeneities of the thermal structure and thus provide im¬ 
portant constraints on the CO distribution. For clarity a contour 
line for T = 5000 K is drawn (white line). The horizontal vari¬ 
ations appear much smaller in terms of logarithmic gas density 
(panel b) which on the other hand exhibits a decrease by several 
orders of magnitude from bottom to top of the model. 

The structuring in temperature and density has its direct in¬ 
fluence on the CO number density uco.td in panel c which is 
alternatively presented as CO abundance (panel d). CO is lo¬ 
cated mostly within a thin layer in the middle photosphere and 
closely follows changes of the gas temperature (see panel a), 
filling the cool interiors of the reversed granulation pattern. Due 
to the steep density gradient in the atmosphere the CO number 
density decreases rapidly with height. Nevertheless, the relative 
number of CO molecules remains high in the upper layers as 
can be seen from the abundance [CO], which is defined as the 
ratio of CO particles and total number of hydrogen atoms (in¬ 
cluding those bound in molecules) on the commonly used log¬ 
arithmic scale with [H]sl2. Here, [CO] is typically close to a 
value of 8.3 but it is strongly reduced at positions where the gas 
temperature exceeds ~ 5000 K (indicated by the white lines). 
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Fig. 2. Single time step of the two-dimensional model (t — 70390 s). The upper row shows gas temperature T (a) and logarithmic 
gas density logp (b), including the velocity field which is represented by arrows. The length of the arrows corresponds to the 
distance a fluid element travels within 15 s time (same scale as x, z-axes). The resulting CO number density ;?co,td is presented 
in panel c and additionally as abundance [CO]td (d) on the usual logarithmic scale with [H] = 12. The third row presents the 
results of the equilibrium calculation (CE) for the same time step (see text for more details): CO number density /?co,ce (e) 
and CO abundance [CO]ce (f). In the next row (g-h) the corresponding results for the UICE calculation are displayed. To 
enable direct comparisons of the different calculations, the data ranges are the same in the second to fourth row. The chemical 
timescale, which has been derived from the CE simulation, is shown in panel i. The differences between the three cases are 
plotted in terms of logarithmic fractions for the pairs TD/CE, TD/UICE, and CE/UICE in panels j-1, respectively. In the lowest 
row all cells with temperatures above 9000 K appear white since those are excluded in the UICE calculations. The corresponding 
color/greyscale-coding is shown next to each panel. The curves mark the height of average Rosseland optical depth unity (dashed) 
and T = 5000 K (solid). The solid curves clearly outline propagating shock waves in the model chromosphere. Note that for the 
cases CE and UICE only the upper part of the model is displayed whereas the whole height extent is shown in the upper panels 
for the time-dependent simulation. 


This is true for the convection zone below the photosphere but 
also for the propagating shock waves which are a ubiquitous 
phenomenon in the upper layers. For the latter the abundance 
reveals a shift of the correlation between temperature and CO 
density which is otherwise strong in the layers below. [CO] is 
still high at the fronts of the displayed shock waves, but it is 
not restored immediately in the wake. This is clear evidence 
that CO does not react instantaneously to thermal changes but 


rather on a chemical timescale (shown in panel i, see Sects. l5~2l 
and 15.41 for more details) which is longer than the dynamic 
timescale due to the wave propagation. The dynamic behaviour 
and in particular the almost instantaneous reaction of the CO 
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Fig. 3. Quantites as function of logarithmic gas density and 
gas temperature for a sample of 139 snapshots: a) log. rel¬ 
ative number of cells, b) log. averaged CO number den¬ 
sity of TD case (log (hco.td)). c) log. averaged ratio of TD 
and CE cases (log («co,td/hco,ce)X d) log. averaged chemical 
timescale (CE); Contours and colors/grey scales are explained 
in the legend right next to each panel. The thick dashed line 
represents the horizontally and temporally averaged stratifica¬ 
tion of the whole sequence, whereas the thick triple-dot-dashed 
line marks a vertical column in the representative snapshot (see 
Fig. ED at x = 2380 km, crossing a prominent shock wave. The 
circles mark heights of z = 250 km, 500 km, and 750 km, re¬ 
spectively. 


number density to temperature changes in the photosphere can 
be seen in an animation provided as online material 3 . 

The results of the corresponding equilibrium calculations 
CE and UICE are shown in panels e-f, and g-h, respectively. 
The logarithmic ratios of the CO number densities of the dif¬ 
ferent cases can be seen in panels j-1 (j: TD/CE, k: TD/UICE, 

3 Additional material can be found at 
http://www.kis.uni-freiburg.de/~sven/research/co.html. 


1: CE/UICE). Note that for UICE temperatures greater than 
9000 K automatically result in a zero CO number density. 
Those cells are blanked in the plots. 

The CE case is very similar to the time-dependent result 
almost everywhere in the photosphere but differs strongly in 
chromospheric shock waves and their wakes. Although the 
[CO] abundance is strongly reduced directly at the wave crest 
in both cases, the time-dependent calculations still result in 
a higher density at the front of the propating wave and in a 
lower one in the wake. This effect is caused by the essential but 
wrong assumption in CE that chemical equilibrium is reached 
before the gas temperature can change, analogous to the ICE 
assumption. Hence, CO equilibrium number density and lo¬ 
cal gas temperature do correlate strongly here, in contrast to 
the time-dependent simulation which does take into account 
finite chemical timescales in combination with the fast propa¬ 
gation of thermal inhomogeneities. In the latter (more realistic) 
case CO is only gradually dissociated when a hot wave arrives 
and does not form again instantaneously after the passage. The 
equilibrium approach by nature cannot take into account this 
behaviour. 

The case UICE yields a very similar picture but with some¬ 
what smaller CO number density. Similar to CE, the devia¬ 
tions from the time-dependent model are small in the pho¬ 
tosphere but become significant in the chromospheric shock 
waves. Finally, we compare CE and UICE: Deviations are very 
small in the photosphere but still get large at high temperatures 
in the centres of chromospheric shock waves. This remaining 
discrepancy must be attributed to differences in chemical in¬ 
put data for which equilibrium constants are used for UICE but 
the chemical reaction network for CE. In particular the reaction 
rate coefficients must be considered as potential source of er¬ 
rors since they are mostly defined for very limited temperature 
ranges only (see Table[2). At the crests of shock waves, where 
the largest discrepancies are found, the temperatures are high 
and thus exceed the stated range for many reactions. 


5.2. Dependence on gas density and temperature 

For a sample of those 139 snapshots that are available both for 
the TD and the CE calculations different quantities are calcu¬ 
lated as function of gas density and temperature. All grid cells 
in this sample - regardless of their spatial and temporal position 
- are sorted into discrete bins for logarithmic gas density log p 
and gas temperature T. The number of matching cells for each 
bin (Fig. H) is highest close to the average stratification (see 
dashed line) but much smaller for small densities and high tem¬ 
peratures in the upper left part of the plot which represents the 
domain of chromospheric shock waves. Furthermore, a grad¬ 
ual bifurcation can be seen with a cool background close to the 
average stratification and a hot component due to shock waves 
above (see also W04). The upper right part of the distribution 
corresponds to the top of the convection zone. 

Panel b reveals the clear tendency of CO being concen¬ 
trated where the temperatures are low and gas densities are 
high. These conditions are best realised in the cooler regions 
above granule interiors in the low and middle photosphere of 
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the model atmosphere. This causes the absolute number den¬ 
sity to correlate best with a temperature of ~ 4400 K. 

More remarkable is the difference between TD and CE, 
which exhibits a significant pattern in panel c. While the ra¬ 
tio of the corresponding CO number densities is small close 
to the average strafication, large values are the rule in the hot 
shock wave domain. The deviations in the wakes, which are 
due to a finite timescale on which the afore dissociated CO 

builds up again, show up as an “eye” around logp- 10 and 

T ~ 3000 K. The triple-dot-dashed line represents a column 
in the exemplary snapshot at x — 2380 km (see Fig.0 which 
crosses a prominent shock wave. Consequently, this line tra¬ 
verses the “eye” and the shock wave domain in the (logp, T)- 
plot. 

Finally, the chemical timescale (Fig. 01) derived in the CE 
calculations exhibits basically three things: First, the timescale 
grows strongly with decreasing gas density. Second, the skew 
in the shock wave domain indicates that there the timescales 
are much shorter due to a reduced final equilibrium concentra¬ 
tion. And third, at a given gas density and with that at a given 
geometrical height a large range of values is present, making it 
difficult to define a meaningful average (see also Sect. 15.4k 

5.3. Height distribution 

Averaging the CO number density horizontally and in 
time over a time span of 50000 s on the geometric 
height scale of the standard model results in the height 
distribution displayed in Fig. 0 It shows a prominent 
peak of max(< nco(z) > x ,t) — 2.81 10 l 2 cm ~ 3 at a height of 
Z — 197 km. Note that averaging on the Rosseland opti¬ 
cal depth scale leads to a somewhat higher maximum 
of max(< nco(r) >* jf ) = 3.68 10 12 cm ~ 3 at an optical depth 
logr = —1.17 which on average corresponds to a geometrical 
height of z = 217 km. The corresponding average abundance 
increases rapidly in the photosphere until a value of ~ 8.3 is 
reached at the bottom of chromosphere where it stays roughly 
constant. For the exemplary time step, which is also shown in 
Fig- El the horizontally averaged CO distribution is very close 
to the corresponding equilibrium values (CE) but starts to devi¬ 
ate increasingly above the photosphere (see thin solid and dot- 
dashed lines). This result is expected since deviations mostly 
occur in the hot shock waves which become prominent usu¬ 
ally in the low chromosphere and above (see Fig. 0). In the 
photosphere the conditions are right for equilibrium. The same 
behaviour is found for the UICE case which perfectly matches 
the average CE profile for the examplary time step. However, 
the maximum is smaller by a factor of 0.86 in the UICE case 
compared to TD and CE (see Fig.0. 

5.4. Chemical timescales 

Already the fact that the CO number density closely outlines 
the hotter regions in the low and middle photosphere (see 
Fig- 0 indicates that there the chemical timescales are small 
compared to the timescales on which the atmospheric condi¬ 
tions change. Strictly speaking, chemical equilibrium refers to 



z [km] 

Fig. 4. Average CO number density as function of geometri¬ 
cal height for the whole time-dependent simulation sequence 
(thick solid) and for the exemplary time step shown in Fig. 0 
(thin solid). Also the results of the equilibrium calculations 
(CE) are displayed for the whole sequence (thick dot-dashed) 
and for the exemplary time step (thin dot-dashed), and for 
UICE calculation (thin dashed, snapshot only). The dotted line 
represents the average for the whole sequence, too, but aver¬ 
aged on the Rosseland optical depth scale. 

a given thermodynamic state of the atmosphere, for instance 
expressed in terms of gas temperature and particle densities. 
For reaching equilibrium it is thus essential that these quantities 
do stay constant. But under continuously varying conditions the 
chemical equilibrium can only be realised if the corresponding 
timescale is shorter than the dynamic one. Otherwise the con¬ 
cept of chemical timescales looses its meaning. 

Nevertheless, chemical time scales are in general very help¬ 
ful since they allow to estimate the importance of individual 
reactions or reaction groups and, furthermore, are needed to 
decide whether the assumption of instantaneous chemical equi¬ 
librium is valid or time-dependent simulations are necessary. 

In principle the chemical tim e scale of a reaction is 
given by its inverse rate (see, e.g., iHerbst & Klemnerer'll973t 
lAvres & Rabinll996l) . for instance, 

Chem = (&7002 «h) * (8) 

for the two-body CO dissociating reaction #7002 (see Table 0. 
The required equilibrium number densities of the involved 
species can be estimated for individual reacti o ns/rea ction 
pairs in the way show n by IHerbst & KlemnereH dl973h and 
lAvres & Rabinl II 19961) but are difficult to be obtained in the 
present case of a complex chemical reaction network. Already 
the assumption of equilibrium number densities for the in- 
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Fig. 5. Chemical timescales calculated from the average rates 
for the three dissociating reactions #14, #4076, and #7002 
(thick solid) and the result from the time-dependent calculation 
in W04 (thin dot-dashed). 

volved species is questionable since those species undergo 
a chemical evolution themselves. Each reaction, which e.g. 
forms CO, does not only depend on the local gas tempera¬ 
ture and thus on the corresponding rate coefficient k but also 
on the number densities of the required reactants. If now one 
of the involved species is not available in sufficient amount, 
owing to another less efficient reaction producing this species, 
this “bottle-neck” reduces the total effectivity of this forma¬ 
tion channel. A chemical timescale of an individual reaction 
is thus only of limited scope under non-equilibrium conditions 
and rather must be seen in the context of the whole reaction 
network. 

The CE calculations described in Sect. !2.3l do take into ac¬ 
count the chemical evolution of the whole reaction network and 
allow to derive total effective chemical timescales and corre¬ 
sponding CO equilibrium densities. For the exemplary snap¬ 
shot the resulting CE timescales are shown in Fig. 0, while 
the dependence on logarithmic gas density and temperature is 
displayed in Fig. IB These two figures make clear that a high 
temperature results in a low final CO number density and a re¬ 
lated short timescale, whereas low temperatures lead to long 
timescales. Hence, the strong thermal inhomogeneities of the 
model atmosphere result in a very large spread of timescales, 
covering many magnitudes, which makes it hard to define a 
representative average. 

Instead, we present average timescales for the three CO 
dissociating reactions #14, #4076, and #7002 in Fig. H] in or¬ 
der to quantify the relative contribution of these reactions. For 
this the reaction rates are calculated for each grid cell for snap¬ 
shots with a cadence of 100 s. After averaging temporally and 
horizontally for each reaction independently, the inverse aver¬ 
age rates give the chemical timescales displayed in the figure. 
Obviously, the time scale for reaction #14 is much shorter than 
those for the other two. Reaction #7002 is typically two or¬ 
ders of magnitude slower, reaction #4076 even around three. 
Consequently the total chemical timescale, given by the inverse 
sum of the rates of the three destruction reactions, is governed 


by reaction 14 alone. At a height of z = 646 km, where the aver¬ 
age mass density is {p) x ,t = 10 15 g cm -3 , the chemical timescale 
for reaction #14 is only ~ 3 s which corresponds to a gas tem¬ 
perature of ~ 4900 K. This can be compared to 0.7 s for the 
rate coefficients given bv lAvres & Rabin (199(» (see Tabled}. 
However, also for this particular column density a large range 
of temperatures and related timescales is found. Hence, the av¬ 
erage numbers should be considered as rough estimates only. 

The timescale for reaction #7002 deviates to some extent 
from the result in W04. This i s m ainly due to the fact that 
the reaction coefficients bvlAvres & Rab in! (119961) were used 
in W04 instead of those bv lBaulch et al.l dl976i) as done in this 
work. The two parameterisations only lead to similar reaction 
rates close to 5000 K and deviate otherwise. Nevertheless the 
average timescales agree quite well above z ~ 700 km. 

5.5. Importance of individual reactions 

Next to the distribution of CO also the relative contribution of 
individual reactions is of great interest. This information allows 
to reduce the chemical reaction network as much as possible 
by excluding irrelevant reactions and thus increasing computa¬ 
tional performance. One also gets to know which reactions are 
most important. 

In our reaction network three reactions are dissociating 
carbon monoxide: 

#14 : CO + H —> C + OH 

#4076 : CO + M C + 0+ M , 

#7002 : CO + H C + 0+ H . 

As we showed in Sect. l5.4l the species exchange reaction #14 is 
by far the most efficient dissociation channel whereas #7002 
and #4076 only contribute 1O 2 and 1 () 3 in the low photo¬ 
sphere and 10~ 3 and 10 4 above, respectively, in terms of num¬ 
ber of dissociated molecules per time. 

There are five reactions that form CO: 


#67 : C + OH -» CO + H 

#104 : O + CH -» CO + H 

#3707 : C + O -» CO + v 

#4097 : C+O + M -» CO + M 

#7001 : C+O + H -» CO + H 


Again, the species exchange reactions are the most effective 
ones under the conditions of the solar atmosphere. These re¬ 
actions, namely #67 (via OH) and #104 (via CH) have simi¬ 
lar rates, although the formation channel via CH seems to be 
slightly faster for temperatures below ~ 7500 K. The other re¬ 
actions, including the direct radiative association (#3707), are 
totally negligible compared to the CH and OH channels. 

Not only the reaction rates but rather the product of the 
rates and the number densities of the involved species is cru¬ 
cial for the change of the CO number density (see Sect. HB 
The average abundance of CH has a maximum of ~ 4.9 in the 
photosphere and decreases above to only ~ 3.9, while OH is 
on average one to two magnitudes more abundant for all atmo¬ 
spheric heights, ranging from [OH]~ 5.0 at the bottom of the 
photosphere to [OH]~ 7.6 in the low chromosphere. Hence, 
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wavelength [nm] 

4662 4663 4664 4665 4666 4667 



Fig. 6. Spectra near A = 4.7 fim with fundamental vibration rotation lines of carbon monoxide (R branch): Average synthetic 
spectra with CO densities from time-dependent 2D model (thin solid) compared to ATMOS3 data (thick solid). Note that the 
UICE calculation exactly coincides with the time-dependent case. 


the formation channel via OH is much more important than the 
CH channel due to the higher number density of the required 
reactant. 

In order to investigate this matter in more detail additional 
runs 4 with different networks have been carried out. In all cases 
we started from the default network presented in Table Q and 
Fig. m In case A all reactions involving CH are eliminated. 
The relative distribution does not change significantly. The ab¬ 
solute values are on average only reduced by ~ 9 % com¬ 
pared to the results with the full network. The same is true 
for the peak of the horizontally and temporally averaged CO 
number density (see Table 0. For case B also the reactions 
#3707, #4076, #4097, #7001, and #7002 are removed. In the 
remaining network carbon monoxide can only be formed or 
destroyed via OH, i.e. reactions #14 and #67. The resulting CO 
number density closely matches the one of the previous case 
outside small regions with extreme temperatures. In contrast, 
the exclusion of OH (case C) produces on average by a fac¬ 
tor of ~ 3 more CO at most positions. In chromospheric shock 
waves, however, nco is not as strongly reduced as in the case 
with the full network. Carbon monoxide is nevertheless con¬ 
centrated mostly in the cool areas of reversed granulation in 
the middle photosphere. Excluding both OH and CH (case D) 
results in a generally smaller CO density, where the largest dif¬ 
ferences are found in shock waves. Again, like in case C, CO 
is not as strongly reduced as for the full network. Note that the 
reduced network is separated into two isolated parts, namely 
H<->Hi and C<-»CO<-»0. The simulation run without the four 
reactions involving the representative metal (case E) leads basi¬ 
cally to the same results as the standard case. For temperatures 

4 Simulation snapshots for the different cases are provided as online 
material. 


below ~ 6000 K the difference in wco between the simulations 
with and without metal is less than one percent and it increases 
only slightly to a maximum of ~ 2 percent for higher temper¬ 
atures. We thus conclude that the adopted catalytic reactions 
and also the abundance of the representative metal are of minor 
importance here. Finally, the network in case F comprises only 
the two reactions for radiative association (# 3707) and colli- 
sional diss ociat ion (# 7002). The same pair was used in W04 
(see also lAvres & Rabin! 1 9961) but here we use the reaction co¬ 
efficients chosen in this work (see Sect. 0. The absolute CO 
number density is on average roughly one order of magnitude 
lower compared to the standard case. The relative distribution 
stays quite similar, except for the upper layers. There the sim¬ 
ple reaction pair produces a smaller relative CO concentration 
that, however, is not as strongly reduced in shock waves as it is 
in the case of the full network. This is in line with the calcula¬ 
tions in W04 for which the fraction of all carbon atoms bound 
in carbon monoxide was only ~ 10 %. This compares to an 
average fraction of 80 % derived with the full network. 


The ratio of the peak of the average CO number density 
of the different cases and the default network in Tabled can 
be used to summarise the aforementioned results. Catalytic re¬ 
actions are negligible (case A). The CH channel is only of 
secondary importance but should not be omitted (case A). In 
contrast #3707, #4076, #4097, #7001, and #7002, including ra¬ 
diative association and collisional dissociation (case F), can be 
neglected (see case C compared to case B). Most important is 
the channel via OH (case B). Excluding it from the network 
(case C) produces fundamentally different results. 
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5.6. CO spectrum synthesis 


6. Discussion 


A sample of 50 snapshots with constant cadence from the time- 
dependent 2D simulation sequence was used as input for the 
spectrum synthesis code RH by Uitenbroek (see Sect. 12.4k For 
all snapshots a spectrum for the same wavelength range near 
4.7 gm was calculated for disk-centre (ji = 1.0). In Fig. El the 
horizontal and temporal average of these spectra is shown in 
comparison with data from the third version of the Atmospheric 
Trace Molecule Spectroscopy (ATMOS3) atlas (IFarme 
iFarmer et al]ll989t http: //remus. j pi . nasa. gov/atmos). 
The ATMOS3 data were scaled to match the continuum inten¬ 
sity of the averaged synthetic spectrum. 


Since the CO number density of the time-dependent and the 
ICE approach closely match in the line-forming layers, also the 
resulting spectra agree almost perfectly. But both exhibit much 
deeper line cores than seen in the ATMOS3 data, similar to the 
results by llJitenbroek ( 1200 01 see his Fig. 11). Owing to the 
dynamics, which tend to be too strong in 2D simulations, and 
the small model extent, the CO number density, the resulting 
depth of the absorption lines, and also the continuum intensity 
vary significantly. Although there are times and positions that 
can produce a better fit than the average shown here, the chosen 
number of 50 snapshots might still be too small to make up a 
representative sample. This is also implied by the fact that the 
average continuum intensity level only corresponds to a bright¬ 
ness temperature of 5400 K. The sample is obviously not per¬ 
fectly representative but rather is too cold. This would produce 
too much CO and consequently too deep line cores. Lines with 
low excitation potential seem to deviate more from the ATLAS 
data than lines with higher excitation potential. Given that low 
excitation lines are formed further up in the atmosphere, this 
trend could point out too much CO in the higher layers. On the 
other hand, the spectra with IC E assum ption show exactly the 
same behaviour, in contrast to ll litenbroekl ( 2000 1 . 


Table 3. Comparison of simulations with different chemical 
reaction networks (case A-F, see text) and the full network 
(REF). Given are the number of remaining reactions n reac , ab¬ 
solute value (max((«co))) and height (z max ) of the maximum 
of the averaged CO number density, the ratio of the peak 
number density and the maximum value of the standard case 
(6nco = max((nco))/max((nco)refX and the gas temperature 
T mc that correlates best with CO number density. 


case 

^reac 

max«n C o» 
[10 12 cm' 3 ] 

£max 

[km] 

Sn co 

Tmc 

[K] 

REF (all) 

27 

3.46 

171 

1.00 

4400 

A (no CH) 

19 

3.15 

185 

0.91 

4400 

B (no CH +5) 

14 

3.15 

185 

0.91 

4400 

C (no OH) 

16 

10.88 

101 

3.15 

5200 

D (no CH, no OH) 

10 

2.37 

129 

0.69 

4700 

E (no M) 

23 

3.45 

171 

0.99 

4400 

F (#3707,#7002) 

2 

0.27 

157 

0.08 

4600 


Spatial distribution. As stated in Sect. 15.11 we find the 
maximum of the average CO height-distribution to be 
2.8- 10 12 cm ' 3 at z = 197 km. Despite differences in the 
assumed rate coeffic i ents, t his is in line with the work by 
lAsensio Ramos et al.l ( 2003 ) who find, judging from their 
Fig. 2, a peak v alue of ~ 3 • 10 12 c m 3 at z ~ 100 km. 
llJitenbroek 1200 01) does not provide a horizontal average but 
the spatial CO distribution in a two-dimensional cut through 
the adopted model (see his Fig. 9). There, the highest values are 
between ~ 3 10 12 cm ' 3 and ~ 10 13 cm ' 3 in the low and mid¬ 
dle photosphere. Next to these high values, which are mostly 
found above granule interiors, also smaller CO concentrations 
are present above intergranular lanes. The horizontal average 
of the distribution is thus in line with the results of this work, 
concerning the maximum values as well as its s pati al distribu¬ 
tion. Furthermore, the findings also agree with lAvresI ( 1981 ). 
who states T 500 < 1 (L 2 as upper limit for the peak of the CO 
distribution. This optical depth corresponds to a geometrical 
height of z = 365 km in our model and is indeed higher than 
the peak found here. Also the fact that CO is mostly located in 
the cool regions of the reversed granu lation pattern matches the 
observed pattern. llJitenbroek OOod) reported on spectroheli- 
ograms in the cores of strong CO lines near A = 4.7 fin 1 that ex¬ 
hibit dark areas as small as 17, lasting for several minutes^sur- 
round ed by bright network-like rims (see also lllitenbroek etHll 
119941 ). The same is not only seen in the two-dimensional model 
presented here but even more clearly in first results from 3D 
simulations that are currently in progress. 

The finding that the bulk of CO is truly located in the pho¬ 
tosphere also fulfills other indirect constraints like the presence 
of 5 minute oscillations. Such prominent intensity oscillations 
were already discovered by In'ovcs & Hal l il972b ) in the core 
of the 3-2 R14 CO line in the fundamental vibration-rotation 
bands_(AV = 1) near a wavelength of A - 4.7 fim (see also 
lAvres & Brauld 1 1 99(ll) . This line should be formed just in the 
high photosphere iNoves & Hall8ll972al) . where 5 min oscilla¬ 
tions can clearly be detected. 

On the other hand, the differences between the synthetic 
spectra and the observed ATMOS3 data might indicate that 
the model still exhibits too much CO - at least in the upper 
layers. One possible cause might be a too simple treatment of 
the radiative transfer in the chromosphere which is still treated 
frequency-independently (grey) and under the assumption of 
local thermodynamic equilibrium (LTE). Already a small er¬ 
ror in the resulting gas temperature can generate significant 
changes in CO number density. Another and maybe more se¬ 
vere source of too high CO number densities might come from 
incorrect or uncertain rate coefficients. However, in order to 
answer these questions a larger and thus more significant set of 
synthetic spectra is needed. 


Necessary modelling. An important result is that CO is 
mostly concentrated within the cool regions of the reversed 
granulation pattern. Hence, particular care should be taken 
for modelling the corresponding mid-photospheric layers. 
Although the basic features are matched, the atmospheric struc- 
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ture in a small two-dimensional model suffers from the vast 
dynamics which are usually exaggerated in 2D compared to 
3D. Consequently, we strongly recommend to proceed to three- 
dimensional models like, e.g., the one described in W04 which 
has been proven to produce a very realistic reversed granulation 
pattern l lLeenaarts & Wedemever-Bohm 120051) . A 3D model 
provides much better statistics which are crucial for the synthe¬ 
sis of meaningful spectra (see Sect. 15.6I> . We expect many re¬ 
sults to remain qualitatively the same, e.g., the dependence on 
temperature and gas density (Sect. l5^2l . the relative contribu¬ 
tion of individual reactions (Sect. 15.5> . and the differences be¬ 
tween ICE and time-dependent approach. In contrast the aver¬ 
age stratification and fluctuations of temperature and CO num¬ 
ber density may be a bit different in 3D compared to 2D but 
still the qualitative picture should be the same. 

Another point concerns the radiative transfer which is 
treated grey, i.e., frequency-independent, in the simulations 
presented here. For comparison we calculate _a short se- 
quence with a non -grey multi-group scheme (cf. lTndwigll992t 
Wedemever 2003). The differences in gas temperature and in 


CO number density are still small in the middle photosphere 
where the bulk of CO is found. In contrast, there are significant 
differences in the chromosphere. The temperature fluctuations 
are much smaller there in the non-grey case. As a result the CO 
number densities of the grey and the non-grey simulation de¬ 
viate most in the hottest parts of chromospheric shock waves. 
While the multi-group approach produces a more realistic ther¬ 
mal structure of the photosphere, the situation is less clear 
in case of the chromosphere. There, neither the grey nor the 
multi-group LTE approach is appropriate. Rather a frequency- 
dependent non-LTE radiative transfer is required which is un¬ 
fortunately not available yet. Nevertheless, the grey approach 
is preferable since it produces results (e.g., shock peak tem¬ 
peratures) that su rprisingly agree better with other works (e.g. 
Carlsson & Stein 1995; Skartlieni ll 998 > than those obtained 


with the multi-group scheme (see W04 for a detailed discus¬ 
sion). But still the modelling of chromospheric radiative trans¬ 
fer remains a major weakness and thus needs to be improved in 
future. 

For the modelling of the number densities of the chemi¬ 
cal species the ICE assumption seems to be valid in the pho¬ 
tosphere where chemical timescales are mostly much shorter 
than the dynamical one so that advection and non-equilibrium 
effe cts can be neglected the re - a result already mentioned by 
(Asensio Ramos et alJ (1 2003 ). This is also implied b y the agree- 
ment of the CO distribution with the results of llJitenbroekl 
( 12000 ). Hence, no time-dependent calculations seem to be nec¬ 
essary for the photosphere which after all contains by far the 
largest absolute amount of carbon monoxide. In contrast, the 
ICE assumption fails close to shock waves in the chromosphere 
where a time-dependent treatment is mandatory. 

The remaining discrepancies in CO equilibrium densities 
between the calculations based on the reaction network (CE) 
and the ICE approach (UICE) can be attributed to differences in 
chemical input data that are equilibrium constants for UICE but 
rate coefficients in the reaction network. In particular the lim¬ 
ited temperature range for which many rates are defined might 
explain why the largest deviations are found at the hot crests of 


the chromospheric shock waves. There gas temperatures often 
exceed the stated range for many reactions. 

Chemical reactions. A thorough analysis of chemical reac¬ 
tion coefficients and additional simulations with altered chem¬ 
ical networks leads us to the conclusion that OH is the most 
important agent for forming and dissociating carbon monox¬ 
ide whereas the CH branch contributes much less - a result 
already suggested by lAvres & Rabini J19961) . In contrast, the 
mere reaction pair composed of radiative association and col- 
lisional dissociation cannot reproduce the observed CO distri¬ 
bution. Catalytic reactions involving a representative metal are 
negligible even if it is assumed to be as abundant as helium. 

Unlike lAsensio Ramos et all (2003) we did not account for 
nitrogen chemistry in our reaction network. In order to check 
the possible consequences of the exclusion of nitrogen, we di¬ 
rectly compared ICE number densities calculated with the RH 
code for the exemplary time step (see Sect.0 for two different 
sets of molecules: First the complete default set (see Sect. E3, 
second a set that only included those molecules that are also 
part of the reaction network used for our time-dependent sim¬ 
ulations. Hence, both sets differ by the following molecules: 
H 2 + , C 2 , N 2 , Oj, CN, NH, NO, and H 2 0. This includes all 
molecules used for nitrogen chemistry by Asensio Ramos et 
al.. The differences in CO number density between the two 
cases are neglible in the atmosphere, in particular at the heights 
where the bulk of CO is present. In contrast, the reduced set re¬ 
sulted in a CO number density lower by up to five orders of 
magnitude at the hot centres of chromospheric shock waves 
but there only little CO is found anyway. Moreover, the aver¬ 
age CO height distribution of our time-dependent simulation 
is very similar to that of Asensio Ramos and co-workers. We 
thus conclude that nitrogen chemistry may be negligible for the 
formation of CO in the solar atmosphere - at least outside the 
problematic hot centres of chromospheric shock waves. 


Chemical timescales. The rates and the corresponding 
timescales span a range of several orders of magnitude, making 
it difficult to define a meaningful average. The results presented 
here should thus serve as a first-order proxy only. 


7. Conclusions 

The presented 2D radiation chemo-hydrodynamic simulations 
produce a spatial distribution of carbon monoxide that agrees 
well with earlier theoretical work and also with observations. 
The synthetic spectra match reasonably well if one consid¬ 
ers the large fluctuations that are unavoidable in the case of 
such a small-size two-dimensional model. More accurate re¬ 
sults can be expected from a 3D simulation which is currently 
in progress. The presented model proved that CO is mostly lo¬ 
cated in the cool regions of the reversed granulation pattern at 
mid-photospheric heights. Some effort was spent to make the 
applied chemical reaction network as realistic as possible on 
the basis of the available chemical input data whose quality is 
quite poor in some cases. Hydroxide (OH) was found to be the 
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most important ingredient for the formation and destruction of 
carbon monoxide. 

The presented comparison between the time-dependent 
simulation and the corresponding instantaneous chemical equi¬ 
librium (ICE) case shows that significant deviations occur only 
near hot shock waves in the solar chromosphere. Hence, ICE is 
a valid assumption for the photosphere that allows to treat CO 
in a simple way in simulations that do not extend beyond the 
photosphere. In contrast, ICE fails in and close to propagating 
high-temperature regions at chromospheric densities. Rather a 
time-dependent treatment, as presented in this work, is manda¬ 
tory for realistic simulations of the solar chromosphere. 

Last but not least we take the agreement of the results with 
other works as demonstration of the validity of the approxima¬ 
tions and algorithms of the upgraded code CO 5 BOLD in the 
case of CO in the solar atmosphere. It now can be used for a 
large range of applications for other stellar types and also other 
chemical reaction networks. 
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